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Abstract 

A method for a time-dependent search for flaring astrophysical sources which can be potentially detected by large 
neutrino experiments is presented. The method uses a time-clustering algorithm combined with an unbinned like- 
lihood procedure. By including in the likelihood function a signal term which describes the contribution of many 
small clusters of signal-like events, this method provides an effective way for looking for weak neutrino flares over 
different time-scales. The method is sensitive to an overall excess of events distributed over several flares which 
are not individually detectable. For standard cases (one flare) the discovery potential of the method is worse than 
a standard time-dependent point source analysis with unknown duration of the flare by a factor depending on the 
signal-to-background level. However, for flares sufficiently shorter than the total observation period, the method is 
more sensitive than a time-integrated analysis. 
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1. Introduction 

The major aim of neutrino astrophysics is to con- 
tribute to the understanding of the origin of high en- 
ergy cosmic rays. A point-like neutrino signal of cosmic 
origin would be an unambiguous signature of hadronic 
processes, unlike y-rays which can also be created in 
leptonic processes. Neutrino telescopes are ideal in- 
struments to monitor the sky and look for the origin 
of cosmic rays because they can be continuously oper- 
ated. The detection of cosmic neutrinos is however very 
challenging because of their small interaction cross- 
section and because of a large background of atmo- 
spheric neutrinos. Parallel measurements using neutrino 
and electromagnetic observations (the so-called "multi- 
messenger" approach) can increase the chance to dis- 
cover the first neutrino signals by reducing the trial fac- 
tor penalty arising from observation of multiple sky re- 
gions and over different time periods. In a longer term 
perspective, the multi-messenger approach also aims at 
providing a scheme for a phenomenological interpreta- 
tion of the first possible detections. 

The search of occasional flares with a high-energy 
neutrino telescope is motivated by the high variabil- 



ity which characterizes the electromagnetic emission of 
many neutrino candidate sources. Resent results ob- 
tained by the IceCube Collaboration fl| indicate that 
high-energy neutrino telescopes have reached a sensi- 
tivity to neutrino fluxes which is comparable to the ob- 
served high energy gamma-ray fluxes of Blazars in the 
brightest states (e.g. the flares of Markarian 501 in 
1997 d and Markarian 421 in 2000/2001 |H). With 
the assumption that the possibly associated neutrino 
emission would be characterized by a flux enhancement 
comparable to what is observed in gamma-rays in such 
states, neutrino flares could be extracted from the sam- 
ple of neutrino-like events with a reasonable signifi- 
cance. 

These astrophysical neutrinos can be searched for 
in several ways. One of the methods for a neutrino 
point source search is to look for events coming from 
a restricted angular region, which could be identified 
with a known astrophysical object. Finding neutrino 
point sources in the sky means to locate an excess 
of events from a particular direction over the back- 
ground of atmospheric neutrinos and muons. These 
events might present additional features that distinguish 
them from background, for example a different energy 
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spectrum or time structure. For sources which man- 
ifest large time variations in the emitted electromag- 
netic radiation, the signal-to-noise ratio can be increased 
by testing smaller time windows around the flare (a 
time-dependent search). In principle there are two ap- 
proaches to neutrino time-dependent searches: 

• Triggered Flare Search: Looking directly for 
photon-neutrino correlations using specific source 
lightcurves from Multi-WaveLength (MWL) ob- 
servations J2l- 

• Untriggered Flare Search: A generalized search 
(MWL data) for neutrino flares, motivated but 
not associated with MWL observations, which are 
scarce and not available for all sources during com- 
plete periods. In addition, there could be a time lag 
between observed photon flares and the associated 
neutrino flares. In the extreme cases high energy 
photons could be entirely absorbed during periods 
of the highest photon production in the source fl. 
This approach however entails a higher trial factor 
penalty than triggered flare searches. As a merit, 
neutrino flares which are not accompanied by ob- 
served electromagnetic counterparts are not auto- 
matically excluded. This approach is also less de- 
pendent on models for the correlation between the 
neutrino and the electromagnetic emission and not 
dependent on the availability of multi-wavelength 
information. 

Here we develop a method that is well suited for an 
intermediate approach in which "periods of interest" 
can be a priori selected on the basis of multiwavelength 
data. The neutrino sample can then be scanned looking 
for significant structures, in a way which is less depen- 
dent on models predicting different correlations with a 
given wavelength. 

An untriggered unbinned flare search was first devel- 
oped and applied to IceCube data, using a compact list 
of pre-defined source directions |Q]|. A time-clustering 
algorithm [1,6], and the unbinned maximum likelihood 
method Oj, |8[] are the basis of this analysis. Such a 
method finds the most significant flare in a long pe- 
riod. The number of trials coming from all combina- 
tions of event times is increased, reducing the signifi- 
cance. However, for flares sufficiently shorter than the 
total observation period, the time clustering algorithm is 
more sensitive than a time-integrated analysis. 

In this paper, we propose an extension of the method 
described in [1]. By including in the likelihood a sig- 
nal term which describes the contribution of many small 
clusters of signal-like events, our algorithm can extract 
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Figure 1: The basic idea of the time clustering procedure. The 
signal-like events are extracted out of a data sample and then 
sorted in time. Each combination of such event pairs defines a 
possible flare time window, which is then tested with a maxi- 
mum likelihood method. 



not only the most significant flare, but also less sig- 
nificant clusters of events distributed over weak flares. 
These weaker flares could be separated by any distance 
in time and will be very difficult to detect or even can 
not be detected with other methods J3, 01 • 

The paper is structured as follows. The algorithm 
is described in Section 2. In Section 3 we apply the 
method to a simulated neutrino search for one flare and 
a few weak flares separated in time. A short discussion 
about the algorithm performance is presented in Section 
4. Conclusions are given in Section 5. 



2. The algorithm 

In this section, we first describe the time-clustering 
algorithm, then we shortly recall the unbinned maxi- 
mum likelihood method [7] and finally we propose our 
new algorithm. 

2.1. The time-clustering algorithm 

The time clustering algorithm [ 1 ] selects the most sig- 
nificant cluster of events in time and returns the mean 
time and width of the corresponding flare. The basic 
idea is shown in Figure Q] In a first step, the method 
selects the most promising flare candidates over differ- 
ent time windows (Atj), which are given by the combi- 
nation of the times of signal-like events from the ana- 
lyzed data set. A signal-like event is defined as having 
Sj/Bj > 1, where 5,- and B, is the background and the 
signal Probability Density Function (PDF) as defined 
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for the time-integrated method [7]2J. Each combinations 
of these event times defines the start and end time (f™ n 
and f™ ax ) of a candidate flare time window (Af ,). For 
each Atj, a significance parameter (the test statistic) TS y 
is then calculated as defined in [7]. Larger values of TS ; 
correspond to data less compatible with the null hypoth- 
esis (i.e. zero expected signal events in the data sam- 
ple tested). Finally, the algorithm returns the best TSj, 
TS max , corresponding to the most significant time win- 
dow over the entire data period analyzed. 

2.2. The unbinned maximum likelihood method 

The unbinned maximum likelihood method J3l de- 
fines the test statistic parameter by: 
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where x s is the source location, h s and y s are the best es- 
timates of the number of signal events and source spec- 
tral index, respectively, which are found by maximizing 
the likelihood £,: 

i*».y.Ao)=n(£Mi-£)4 (2) 

/=i 

where is the number of all events in the data sample. 
The background PDF is given by: 
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where f>^ pace describes the space distribution of events at 
a given region of the sky with the solid angle dQ., p ertelgy 
is the energy distribution and p Ume describes the back- 
ground time distribution. These PDFs can be calculated 
purely from data. In general, due to analysis cuts, Earth 
absorption effects and the detector geometry, the spa- 
tial probability p s . pace and the energy probability f™ elgy 
depend on zenith, f?,, and azimuth, </>,. 
The time probability P ame is defined by: 
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where AT fata is a normalization constant for the whole 
data taking period. 

The properties of signal events are taken from a dedi- 
cated Monte Carlo (MC) signal simulations. The signal 
PDF, Si is given by: 

S t = P; pace (| 2 i -2,\,o-dI? aa (E i ,6 i .Y,)rf a> , (5) 



1 To calculate the ratio, 5,/B, > 1, only the spatial and energy 
terms in the PDF's are included. 



where the spatial probability p spdce is a Gaussian func- 
tion of | x, ■ - x s |, the space angular difference be- 
tween the source location x s and each event's recon- 
structed direction x/, and <r, the angular error estimate of 
the reconstructed track. The energy probability p energy ? 
constructed from signal simulation, is a function of the 
event energy estimator the zenith coordinate f?,, and 
the assumed energy spectral index of the source y s (such 
that % ~ E~ 7s ). Pf me is the time probability which 
generally is a constant value (i.e. taken to be uniform 
in time) if no flare structure is assumed. For each time 
window tested (Af/ = f™ ax - f™ n ), the time probability 
is given by: 
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where f, is the arrival time of the i* event and H is the 
Heavyside step function. Note, that by using this def- 
inition for the time probability in the likelihood X., we 
count only events which fall inside the current time win- 
dow Atj (i.e. the signal PDF S ,■ is zero outside the se- 
lected time window). 

As shown in ||8|] the unbinned likelihood method will 
preferentially find shorter flares, making it less power- 
ful for flares of durations longer than roughly one day. 
The solution to this problem is to use a marginalization 
term (^) in the likelihood £ [8]. This gives a more 
uniform exposure to find flares of different widths and 
leads to a redefinition of the test statistic: 
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2.3. A method to search for multiple flares 

We propose here some extensions to the above men- 
tioned procedures to identify a series of weak flares by 
incorporating this information into the likelihood func- 
tion. This is done by restricting the search to dou- 
blets of signal-like events, and using the value of the 
test-statistic of these individual flare candidates as their 
weights in a stacking-like calculation of the global max- 
imum likelihood. 

First, we extract all doublets that can be formed out 
of all signal-like events (S j/Bj > 1) over the entire data 
taking period A7d ata . This step serves to isolate all pos- 
sible (and smallest) time windows that compose the sig- 
nal contribution in the tested data sample (the total num- 
ber being M). We call these time windows "data seg- 
ments". Note, that by construction, there could be a cer- 
tain degree of overlap between different data segments 
if larger event multiplicity will be studied, see Figure |2] 
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Top. The choice of doublets is physics-motivated, be- 
cause for a neutrino detector like IceCube, the signal ex- 
pectation is not much more than a few signal neutrinos 
per year from the strongest astrophysical sources [0. 
Moreover, we focus on weak multiple flares. 

Then for each time window Atj a minimization of 
\og(X.) as defined in Eq. (2) is performed, with n s and 
y s as free parameters, and the individual test statistic 
TSjUr is calculated. This step serves to estimate the 
possible signal contribution in each data segment. All 
time windows are then sorted according to TSjIa;,, as it 
is shown in Figure [2] (Middle). Some of these data seg- 
ments will contain real signal events and some of them 
are likely due to background fluctuations. Our aim is 
to extract the optimal (best suited) number of data seg- 
ments (M opt ) which compose the total signal contribu- 
tion injected in the overall period A7d ata . 

For this purpose, we propose a modification of the 
single-source likelihood function (Eq.(2)) by including 
a signal term (Sf'im)) describing the contribution of 
each data segments: 

N 

l{n s ,y,m) -\\i^Sf\m) + (l - (8) 

where m is the number of data segments in to which the 
overall signal contribution can be decomposed (a free 
parameter m < M) and N is the total number of events 
in the time period Atta- 
in other words, in order to include the contribution 
to the signal from multiple flares, the one-source signal 
term S , is being replaced by the sum of signal sub-terms 
over m data-segments: 



Z7 = i WixSjQjti-tslEi^Atj) 



(9) 



where W ! is a weight which describes the strength (sig- 
nificance) of the doublet contained in each data seg- 
ment. This weight should be proportional to the ex- 
pected number of neutrino events seen in the detector. 
As we will show later the test statistic is quite well cor- 
related with the true number of injected signal events. 
Thus we take W j = TS 7 -| A « r 

The likelihood given in Eq. ([8]l in combination with 
Eq. (9) is usually used in stacking searches for differ- 
ent types of astrophysical point sources ifioll . but here 
the individual sources are data segments of the signal- 
like events found in the data from the same object. 
The stacking method uses the fact that signal and back- 
ground increase differently when observations for mul- 
tiple sources (or flares) are being added. Stacking multi- 
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Figure 2: Sketch of the flare stacking procedure. Top: A data 
set is divided into segments being defined as in FigurefJJ Some 
of these segment will include signal contribution, others back- 
ground fluctuations only; Middle: Sorted values of the test 
statistic for each data segment from a maximization of the 
likelihood in Eq. (2) as a function of the data segment index 
m. Bottom: Evolution of the test statistic from a maximiza- 
tion of the likelihood in Eq. (8) as a function of the number 
of data segments being stacked following Eq. (9). The maxi- 
mum value of this parameters defines the optimal number of 
data segments to be stacked M opt . 



pie flares can therefore suppress the overall background 
and hence can increase the signal sensitivity. 

In order to estimate the optimal number of data seg- 
ments M opt for a given number of m segments (starting 
from m — 1) we minimize the log( J £(n s , J, m)) with n s 
and y s as free parameters. The minimization returns the 
best estimates for the number of signal events h s and for 
the spectral index of the source y s , and the "global" test 
statistic is calculated from: 



TS(m) = -2 log 



X.(x s , n s = 0) 



£(x s ,h s ,y s ,m) 



(10) 



Then, the optimal number of data segments to be 
stacked (M opt ) is chosen according to the maximum of 
TS(m) (see Figure|2]Bottom). As a result, the following 
important parameters are extracted: 

• M opt : the optimal number of segments which com- 
pose the signal contribution over the time AFdata- 
The optimal set could be either subsequent in time, 
forming one significant cluster of events (one flare) 
for a given source location, or theses segments can 
be distributed among a few (sometimes less signif- 
icant) flares separated in time. This is shown in 
Figure|2](Top). 
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• h s : the total number of expected signal events 
summed over the M opt individual segments. 

• y s : the spectral index of the source (the flare), as- 
sumed to be the same for each M opt segments. 

• TS(M opt ): the maximum value of the test statistic 
calculated with the modified likelihood function of 
Eq. QOll. 

• Ar(M opt ): the flare duration calculated for M opt 
segments. This is here defined as the time between 
the start time of the first data segment and the end 
time of the last data segment. 

The overall significance of the optimal configuration 
M opt can be determined using MC simulations by apply- 
ing the same procedure to a large number of simulated 
data samples. This will automatically account for ef- 
fects of trial factors. The trial factors arise from testing 
different time windows for the same source direction. In 
order to represent a background-only observation, the 
properties of the data events (e.g. zenith, azimuth, time, 
reconstruction error and energy estimator) are sampled 
from their distributions. Data of neutrino experiments 
with a signal flux can be simulated by injecting signal 
events on top of the background data. 

3. Results 

The method described in Section 2 has been ap- 
plied to 10,000 MC background samples (scrambled sky 
maps). Background events were randomly chosen re- 
producing realistic PDF's (energy, angular resolution) 
extracted from [9]. More precisely, the number of back- 
ground events in a given declination band (n^ nd ) is es- 
timated from «^ nd = n bg X ^ff X f^ff , where « bg is 
the number of background events in a bin and A Da nd and 
Abi n are the solid angles of the band and bin, respec- 
tively. Assuming i% g = 2.5 background events for decli- 
nation 6 = 15° as in [9] and AT^ dta = 40 days we are left 
with 355 events distributed in a declination band of size 
±6°. The angular reconstruction error cr, of each event 
was generated based on the cumulative spread function 
from [9] with a median resolution of about 0.8° and as- 
suming that o-j does not depend on the energy of the 
reconstructed event. This assumption in general may 
not be true, but it won't affect the results of this test. 
The same number of scrambled maps were generated 
for signal events injected on top of the background (MC 
background plus signal simulations). A neutrino point 
source at declination 15° was considered. Signal events 
were injected around this declination d s along a band 
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Figure 3: Sketch of the three examples discussed in this work 
i.e. one flare (Example 1) two individual flares (Example 2) 
and three weak flares separated in time (Example 3). 

6 S - 6° < 6 S < 6 S + 6° o an d with azimuth randomly cho- 
sen from 0° to 360°. Individual event directions were 
smeared according to the Point Spread Function (PSF) 
from [9] with a median angular resolution of about 0.8°. 
We simulate signal strengths following Poisson statis- 
tics with mean values of 4, 8, 9 and 12 signal events. 
The energy PDF follows what was found for an E~ 2 
energy spectrum in H. For each MC signal sample 
(scrambled sky map), signal events were randomly in- 
jected inside a time window defined by various flare du- 
rations, in the range from 0.01 day to 30 days. The most 
of the results presented in this paper was obtained con- 
sidering a total period of Ar^ta = 40 days. This is mo- 
tivated by the fact that an obvious application of this 
method is to test periods of interests, preselected with 
the help of multi-wavelengths data. 

We tested the method for the three following cases 
which are also illustrated in Figure [3] 

• Example 1: Signal events (with Poisson-mean=8) 
are injected inside a 9 days time window. The start- 
ing point of the flare is randomly chosen over a 
data taking period A7d a ta of 40 days. An exam- 
ple of a physics case corresponding to this config- 
uration is the reported y-ray flare from the Crab 
Nebula HI], iQj], 0]. Note, that this case would 
correspond to a rather "strong" flare. As we can 
see for example from [8], the number of events re- 
quired for the discovery of such a flare is larger 



2 Events outside this band hardly contribute to the signal PDF since 
they are too far away from the source compared to the point spread 
function of (gj- 
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Figure 4: (A), (B) Distribution of the best fit values for the number of signal events h s and the source spectral index y, and their 
correlations for background-only simulations (Left) and with 8 (Poisson-mean) signal events injected on the top of the background 
(Right) for a source at declination 15°. (C), (D) The optimal number of data segments M opt as a function of the duration of the flare 
AT (i.e. the maximum time difference between any data segments yielding the best fit configuration M opt ). 



than 7 in case of a search with an unknown burst 
duration. 

• Example 2: The total number of signal events 
(Poisson-mean=8) is the same as in Example 1, 
but individual events are injected over two time 
windows of duration 4.5 and 9 days, respectively 
and 22 days apart. The average number of injected 
signal events is the same for each flare (Poisson- 
mean=4). The expected number of events required 
for the discovery depends on the flare duration and 
it is smaller for the flare with shortest duration. In 
other words we can think of this example as a case 
with a "strong" flare and one weaker flare. These 
two flares are separated in time, but they can also 
be treated as one long period of enhanced emission 
with a duration about of 35 days. Standard point 
source search methods applied to this case will find 
only the strongest flare i.e. the first one. This is 
because standard methods use the one-source like- 
lihood function given by Eq. (j2]l and can look only 
for a minimum corresponding to the most signifi- 



cant cluster of events from a given source location. 

• Example 3: The total number of signal events 
(Poisson mean=8) is the same as in previous exam- 
ples, but individual events are injected over three 
time windows with duration of 4.5 days, 4.5 days 
and 9 days, respectively. Each flare is simulated 
with a similar strength with Poisson-mean 3, 3 
and 2, respectively. This example describes three 
"weak" flares. Note that the number of events re- 
quired to claim discovery for each of them at a So- 
level is larger than 4, as we can see from 0. In 



3 In this work we injected signal events according to a uniform dis- 
tribution while in 1 8] a Gaussian distribution was considered. While 
a general comparison is still possible in some cases the exact numeri- 
cal results are different. The duration of a time window, for example, 
is somewhat different. For a uniform distribution with non zero val- 
ues in the interval At the standard deviation is defined according to 
cr = At/ Vl2 with <x denoting the standard deviation of the Gaussian 
distribution. Thus 4.5 and 9 days width of the flare corresponds to 
(T; — c about 1.3 and 2.6 days or to 0.0035 and 0.0070 year, respec- 
tively. Thus for the flare with a Gaussian mean time <t, = 0.0035 we 
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other words such individual weak flares cannot be 
found at a 5<x level by the standard point sources 
search methods. 

3.1. Example 1: one flare 

In Figure 3] the distribution of the most important pa- 
rameters i.e. the number of signal events h s , the spec- 
tral index y s and the optimal number of data segments 
(i.e. signal-like doublets) M opt and their correlations 
as obtained by MC simulations are presented for the 
background-only simulations and the background plus 
signal cases. The number of injected signal events fol- 
lows a Poisson distribution with a mean of 8. Events 
are injected over a time window of 9 days duration. The 
pure background case illustrates how signal-like flares 
can be mimicked by pairs of events (sometimes also 
triplets) distributed over very short time windows with 
durations of less than a day. The signal plus background 
plots show that the proposed method finds the right flare 
i.e. it recovers the true values of the spectral index (2), 
the number of injected signal events (Poisson-mean 8 
events) and the flare duration^ (9 days). In addition, the 
flare can be decomposed into about 8 data segments, see 
Figure|4](D). 

In Figure|5](A) the correlation of the global test statis- 
tic TS(M opt ) as a function of the number of injected sig- 
nal events n,(MC) is shown. The solid line represents 
a linear fit to the mean values of the test statistic cal- 
culated for a fixed number of injected signal events. A 
linear correlation is found between the test statistic and 
number of injected signal events. This justifies using the 
test statistic as a weight in the modified likelihood £. of 
Eq. (9). We also checked the results obtained when all 
weights are fixed to one as well as to the square of the 
test statistic. We find that such modifications lead to 
slightly worse results i.e. we observed a 5-10% worse 
agreement between the true and the estimated values of 
h s and y s . 

In Figure[5](B) the number of signal events is shown 
as a function of the optimal number (M opt ) of data seg- 
ments which compose the flare. As expected, the num- 
ber of signal events increases when more segments are 
being added. The distribution shows a maximum cor- 
responding to the true number of injected signal events 
(the true Poisson-mean of 8). The inset in Figure[5](B) 
shows the Gaussian fit to the estimated number of sig- 
nal events, h s with a mean value of 8.2, which is good 
agreement with the true value at a level of about 2.5%. 



need about 4 events for discovery using the method labeled by "As- 
sumed Burst Time (Energy)" in Figure 4 of |8[. 

4 The exact duration of the flare can be slightly less since the signal 
events were randomly injected inside the 9 days time window. 



The proposed algorithm can therefore recover the 
most important parameters characterizing flares, like the 
duration, the source spectral index and the number of in- 
jected signal events with uncertainties not larger than a 
few percent. In other words the proposed algorithm can 
"decompose"' the flare into small data segments which 
contain signal-like doublets and can effectively reject 
background events from the entire data period. 

3.2. Example 2: two flares separated in time 

In Figure [6] (A) we show the distribution of the op- 
timal number of data segments M opt as a function of 
the corresponding flare duration AT as defined in sec- 
tion 2.3. In Figure [6] (B) the best fit spectral index y s 
is shown as a function of the estimated number of sig- 
nal events h s . Also in this case the proposed algorithm 
recovers the overall flare duration, the source spectral 
index and the number of injected signal events. How- 
ever, for this particular case, the differences between 
the true and the estimated values of the physics param- 
eters obtained from the minimization are slightly worse 
(about 5%). This is because for the larger flare durations 
we should expect a larger contamination of background 
events. 

In Figure [6] (C) the distribution of the mean time 
Tqj = tj+Atj/2 calculated for each (/ h ) signal-like dou- 
blet is shown. An accumulation of signal-like doublets 
for a time period below 5 days and between 25 and 35 
days is visible, corresponding to the first and the second 
flare, respectively. The proposed algorithm can there- 
fore find not only the most significant flare, but also the 
weaker one separated in time. Such a flare would not be 
detected by other existing methods (like fl). 

By repeating the minimization with the modified like- 
lihood £. for a fixed number of m data segments sorted 
in time (starting from m = 1 and ending with m = M opt ), 
we can get more information about the "internal" struc- 
ture of the flare. For example, in Figure 0(D) changes 
of the global test statistic as a function of the flare dura- 
tion AT(m) calculated for m < M opt data segments is 
shown. The test statistic increases when more segments 
are added in time, but finally reaches a saturation for 
the time corresponding to the overall period of enhance- 
ment emission (35 days). This behavior is strongly cor- 
related with the distribution presented in Figure |6](C). 
For the first flare (below 5 days) and second flare (from 
25 to 35 days) the increase of the global test statistic in 



5 For a fixed number m of data segments, the flare duration is de- 
fined as the time between the start time of the first in time data segment 
and the end time of the last data segment. 



7 




n s (MC) n s 

Figure 5: (A) Correlation of the test statistic and the true number of injected signal events. The solid line represents a linear fit to 
data. Full circles correspond to an average value of the test statistic for a fixed number of injected signal events. (B) The optimal 
number of doublets which compose the flare with 9 days duration as a function of the average number of expected signal events. 
The inset shows the distribution of the best fit values of the number of signal events which can be well fitted with a Gaussian. 



time is much larger than during the period correspond- 
ing to the 22 days time gap. During this last period 
larger variations can also be observed, which is an indi- 
cation for background fluctuations. A similar behavior 
is seen for the average number of signal events (h s (m)) 
as a function of the flare duration (Figure [6] (E)). Note 
that both distributions presented in Figure 0(D) and (E) 
are averaged over different MC realizations. 

For completeness, Figure [6] (F) shows the average 
spectral index (y s (m)) for a large number of different 
MC realizations as a function of time. The average best 
fit spectral index is about 2.2 and differs by about 8.5% 
from the true value (y s = 2). Similarly to the previ- 
ous distributions, the fluctuations of the spectral index 
are strongly reduced during the periods corresponding 
to the injected signal. The effect is particularly visible 
for the second flare, between 25 days up to 35 days. 

We conclude that the proposed algorithm can recover 
the global parameters describing the flares and pro- 
vide additional information about the time development, 
their duration and the number of signal events per indi- 
vidual flare, even when they alone are below the thresh- 
old for detection. 

3.3. Example 3: three weak flares separated in time 

In this case we simulated three individual weak flares. 
Each flare cannot be individually detected at a 5<x level 
by using the standard point source search methods. In 
principle, if all events could be injected into a single 
flare of duration about 28 days (cr, - 0.02 year) the al- 
gorithm proposed in [8] would yield a discovery. For a 
flare of such a duration, more than about 7 events are 
needed for discovery at a 5<x level using the method 



labeled by Assumed Burst Time (Energy) in Figure 4 
of fl. However, the flare will be found only if these 
7 events will form one cluster of events compact in 
time. This is because the unbinned maximum likeli- 
hood method with a single-source likelihood function 
can only search for a maximum which corresponds to 
the most significant flare (cluster of events) from a point 
source. This is certainly not the case for this example, 
where 8 events are distributed among three individual 
and well separated flares. Thus, such a structure cannot 
be detected at 5cr level by the standard method. 

Figure [7] (A) and (B) show that also in this case the 
method can recover the true values of the spectral index 
(y = 2), the mean number of injected signal events (8) 
and the total flare duration (28 days). In addition, the 
overall signal injection can be decomposed into about 8 
data segments (doublets). 

Sorting these data segments in time we obtain the dis- 
tribution of individual flares in time as shown in Fig- 
ure[7j(C). Three flares can be clearly distinguished, with 
a duration of about 5 days, 5 days and 10 days, respec- 
tively. 

The average value of the test statistic and the average 
number of signal events from the best fit (of the likeli- 
hood £) is presented in Figure [7J(D) and (E). We ob- 
serve a similar behavior as in the Example 2. However, 
due to a smaller time gap between individual flares, a 
factor 4 smaller than for the flares studied in the Exam- 
ple 2, the structures are less pronounced. Note that the 
distributions presented in Figure |7](D) and (E) start to 
saturate at a point corresponding to the overall period of 
enhancement emission (above 28 days). 

Finally, in Figure [7j (F) the average spectral index 
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Figure 6: Example 2: Results of the search for two individual flares separated in time using the method presented in this work. (A) 
Distribution of the optimal number of doublets (M opt ) as a function of the total flare duration AT; (B) The spectral index versus the 
expected number of signal events; (C) The mean time of each data segment; (D), (E), and (F) The distribution of (TS(m)j, (h s (m)} 
and spectral index (y s (ni)) as a function of the flare duration AT(m). Also shown are l-<x ranges. The vertical dashed line indicates 
the overall period of enhanced emission (35 days). See the text for more details. 



{y s (m)) as a function of time is shown. The source spec- 
tral index has a value of about 2.2 and differs by about 
10% from the true value (y s = 2). Fluctuations of the 
spectral index also increase for times greater than the 
total flare duration (above 28 days). 

This example demonstrates that the proposed algo- 
rithm can find a few weak flares separated in time, 
which cannot be found by other standard methods . 

4. Performance of the method 

In this section we calculate the number of events 
needed for discovery (i.e. the discovery potential) for 
the cases of Example 1, 2 and 3 and for different flares 



durations AT and overall data periods A^ata- The dis- 
covery potential is defined as the average number of 
signal events required to achieve a p-value less than 
2.87 x 1(T 7 (one-sided 5<x) in 50% of the trials. The 
comparison below gives us an idea about the perfor- 
mance of the proposed method and its limitations. 

In Figure [8] we show the distribution of the global 
test statistic TS(M opt ) for MC background-only simu- 
lations (A) and for MC background-plus-signal simula- 
tions (B), with signal events distributed according to a 
Poisson distribution with mean 4, 8 and 12 (the case of 
Example 1: one flare with duration 9 days). From Fig- 
ure[8](A) we can estimate the TS|5o-(M opt ) threshold for 
a 5cr execss in the background only case, which is 101. 
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Figure 7: Example 3: Results of the search for three individual flares separated in time. (A) Distribution of the optimal number 
of data segments (M opt ) as a function of the total flare duration AT; (B) The spectral index versus the expected number of signal 
events; (C) The mean time for each data segments; (D), (E), and (F) The distribution of (TS(m)j, (h s {m)) and the spectral index 
(y s (m)), respectively. The vertical dashed line indicates the overall period of enhanced emission (28 days). See the text for more 
details. 



The 5<x threshold was calculated for a band of ±6°Hcen- 
tered at declination 15° and for an overall data period of 
AFdata = 40 days. From Figure [8] (B) we can see that 
this threshold is passed in 50% of trials when the aver- 
age number of added signal events is 8. The discovery 
potential is therefore about 8 events. 

In Figure [9] (A) the global test statistic distribution 
for the case of Example 1, 2 and 3 is shown. The me- 
dian is 101, 96 and 91, respectively, which means that 
the number of events required for discovery (assuming 
a threshold of 101) is 8.0, 8.4 and 8.9. Note that we 
have a similar number of events required for discovery 



Calculations with a larger band ( ±10°) does not show changes in 
the 5<x threshold. 



for Example 2 and Example 3 in which signal events 
are injected over a few flares, but such cases cannot be 
tested with the method presented in fH. 

In Figure |9](B) the discovery potential as a function 
of different overall data periods is shown. The number 
of events required for discovery as expected depends on 
the overall data period A7d a ta considered. This is be- 
cause the number of signal-like events increases with 
time and therefore also leads to a higher 5<x threshold. 
For example, for an overall data period of 40 days the 
average number of signal-like events is 5, while it is 43 
for a data period of 365 days. Such changes in the num- 
ber of signal-like events lead to an increase of the 5cr 
threshold from 101 to 180. As a consequence the num- 
ber of events needed for discovery increases by about 
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Figure 8: (A) Probability distribution of the global test statistics TS(M opt ) for background-only simulations. The solid line shows an 
exponential fit to the distribution used to obtain the 5<x threshold for the background only case; (B) Distributions of the test statistic 
TS(M opl ) for background plus signal simulations, where signal events are injected according to a Poisson distribution with mean: 
4, 8, 12 signal events, respectively. The vertical line corresponds to the 5<x threshold. 




Figure 9: (A) Distribution of the test statistic TS(M opt ) for background plus signal simulations with a Poisson mean of 8 signal 
events. The results obtained for one flare with duration 9 days (Example 1), two flares with a total duration (including gaps) of 
35 days (Example 2) and three weak flares with a total duration of 28 days (Example 3) are compared. The vertical dashed line 
indicates the 5cr threshold for background only case; (B) Discovery potential as a function of an overall data period Ar d!lt!1 for one 
flare with different durations and for two flares and three weak flares. The dashed line shows results of a time-integrated analysis 
performed in this work with an energy term in the PDF. 



63%, from 8.0 up to 13 for a flare with a duration of 9 
days. Comparing with the results presented in [8] where 
we need about 9 events for discovery for a method us- 
ing unknown burst time with energy PDF and one year 
data period (Figure 4 in H), we see that our method re- 
quires about 44% more events for discovery in case all 
events are injected in one single flare (Example 1). This 
is because our algorithm stacks also background fluctu- 
ations, and thus leads to a higher 5<x threshold than the 
threshold obtained by a single-source likelihood based 
method. However, it must be noticed that this compar- 
ison is only qualitative, since the signal-to-background 
ratio, which reflects the detector efficiency and signal 
properties used in both simulations, is different in the 
two cases. A more quantitative comparison of both two 



methods has to be achieved on the same simulated data 
set. As we also expect, the number of events needed 
for discovery decreases, when we consider flares with 
shorter duration. For example, for a flare with dura- 
tion of 30 days, 1 day, 0.1 day and 0.01 day in aver- 
age about 9, 6, 5 and 4 events are needed for discovery 
within A7d ata = 40 days. Note that in this case for flares 
with shorter duration (below 0.1 days) the number of 
events is smaller than the number of events for a time- 
integrated analysis, see Figure|9](B) for the method la- 
beled "Time Integrated Analysis (Energy+Space)". 

In Figure [9] (B) the results of calculations for condi- 
tions of Example 2 and 3 are also shown. In general, we 
see a similar trend i.e. the number of events needed for 
discovery increases when a larger overall data period is 
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studied. It is also seen that for Example 2 the number 
of events for discovery is smaller than for a single flare 
with similar duration. However, care must be taken in 
such a comparison because, as we already pointed out 
before, multiple combinations of individual flares can- 
not be found by the standard point source search algo- 
rithm as proposed in [8]. This is especially true for Ex- 
ample 3. In other words, the performance of the method 
improves when we study more flares distributed in time, 
especially if they are weak. 

Up to now we also do not exploit the fact that the 
algorithm provides us information about the number of 
individual flares distributed in time or even about the 
number of data segments which compose one flare or 
a few individual flares separated in time. Using such 
information, we can calculate a better sensitivity. By di- 
viding the sensitivity by the number of flares or data seg- 
ments M opt , we can calculate a sensitivity per individual 
flare or even per segment. For examples presented here 
(which correspond to signal-like doublets) the number 
of M opt is about 8 and the number of individual flares is 
larger than 1 , hence the sensitivity can be improved by 
about the same factor. 

5. Summary & Conclusion 

We have presented a method to search for neutrino 
flares from point sources without an a priori assumed 
time structure. The method considers only data seg- 
ments which contain signal-like doublets, and uses a 
test-statistic term as their weights in a stacking-like cal- 
culation for the global maximum likelihood. We have 
shown that this method can recover the true values of 
the source spectral index, the flare duration, and the to- 
tal number of injected signal events within uncertainties 
not larger than 10%. In addition, our algorithm pro- 
vides relevant physics information about the distribu- 
tion of flares in time and their internal structures. This 
information can be used to calculate a sensitivity per in- 
dividual flare, which is usually better than the sensitivity 
obtained from the other methods. 

For standard cases (one "strong" flare), the discov- 
ery potential of the method is about 44% worse than 
the standard point source analysis with unknown dura- 
tion of the flare. This is because our method stacks all 
significant background fluctuations in a given period of 
data and leads to a higher threshold for discovery. How- 
ever even in such a case, the number of events required 
for discovery is smaller compared to a time-integrated 
search as soon as the flare duration is less than a few 
hours. When the number of individual flares analyzed 
is increased the number of events needed for discovery 



decreases, especially in the case of a few weak flares 
distributed over a longer period (from a few days up to 
100 days). Such cases of a few weak flares cannot be 
discovered by the standard point search algorithm [H. 
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